# ******** Codes for the paper "Distortionary Effects of Direct Lending on Firm Quality Dynamics" *******
# Note: this file will be called by "MasterFile.R". 
# remove(list=ls())
# setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

# ------------------ Plot Specifications ---------------------
margins = c(3.5, 3.5, 1, 1)
plot_slippery_slope_two_economy = TRUE
plot_width = 4.2
plot_height = 3.8

# # ------------------ Set of initial parameters -------------
# Initialize parameter values 
par_init = data.frame(AH=0.55, theta=9, beta=0.4, lambda_F=6, avg_zeta=0.15, d_bar=0.25, g_bar=0.14, 
                      r=0.066, delta=0.2, lambda=0.06, eta_H=0, eta_L=0.062, w0=0.5, AH_AL_ratio = 3.7)
par = par_init
# calibration function 
source("Model scripts.R")  # This script solves the model and provide simulation functions 
result=calibration_moments(par)
equi_sol = result$equi_sol
par = equi_sol$par

# ------- Calibration: solve the nonlinear equation system -----------
moment_targets = c( 0.1, -0.093, 0.49, 0.45, 0.1, 0.36 )
par_index_calibration = c(1:6)
# These moment targets are for (1) investment/capital, (2) avg output drop in crises, (3) creditor_recover, 
# (4) avg productivity to capital ratio,(5) bankruptcy reduction due to govt intervention, (6) private-sector debt/GDP ratio
equi_sol$key_moments
equi_sol$d_bar
solve_moments <- function(par_vec){
  par[par_index_calibration] = par_vec
  calibration_results=calibration_moments(par)
  residual = as.numeric(calibration_results$key_moments - moment_targets)
  print( paste0("residual error = ", round(sum(abs(residual)),digits=6)  ) )
  return(residual)
}
calibration_sol = fsolve(solve_moments,  as.numeric(par[1:6])  )  # Calibration is finished!!!
# Note: we cannot use nleqsv here because it cannot support nested solution. 
par[par_index_calibration] = calibration_sol$x
calibration_results=calibration_moments(par)
equi_sol = calibration_results$equi_sol
moments = cbind(calibration_results$key_moments, calibration_results$other_moments)
print("moment residuals:")
print(calibration_sol$fval)
print("moments:")
print(moments)
print("par values:")
par
Sys.sleep(5)
print("other moment:")
calibration_results$other_moments

# Output the calibrated values
save(par, file = "baseline_calibration.Rdata")

# Then assign values for next processing
g_bar = par$g_bar; d_bar = par$d_bar; dt = equi_sol$dt;
AH = par$AH;  AL = par$AL
lambda=par$lambda; beta=par$beta; r = par$r; eta_H=par$eta_H; eta_L=par$eta_L
w_grid = equi_sol$w_grid
f_list = calibration_results$f_list
solve_individual_optimal_faster = f_list$solve_individual_optimal_faster; F_fun = f_list$F_fun; Delta_output_fun = f_list$Delta_output_fun; solve_equilibrium=f_list$solve_equilibrium; Simulate_model=f_list$Simulate_model; solve_other_moments=f_list$solve_other_moments; 
muK_fun=f_list$muK_fun; DeltaK_fun=f_list$DeltaK_fun; muw_fun=f_list$muw_fun; Deltaw_fun=f_list$Deltaw_fun; investment_cost_fun=f_list$investment_cost_fun;  l_dN_fun=f_list$l_dN_fun;  post_crisis_fraction_of_capital=f_list$post_crisis_fraction_of_capital; output_drop_and_gbar=f_list$output_drop_and_gbar;  
solve_welfare = f_list$solve_welfare
simulation_fun=f_list$simulation_fun
avg_zeta = par$avg_zeta; 
qL = equi_sol$qL; qH = equi_sol$qH;  

equi_sol_baseline = equi_sol    # This is the baseline equilibrium solution 

# Basic moments
equi_sol_no_govt = solve_equilibrium(d_bar, g_bar=0)
output_drop_and_gbar( par$g_bar,  equi_sol$w_avg, qH = equi_sol$qH, qL=equi_sol$qL  )
output_drop_and_gbar( 0,  equi_sol$w_avg, qH=equi_sol_no_govt$qH, qL=equi_sol_no_govt$qL  )
Deltaw_fun( equi_sol$kappaH, equi_sol$kappaL, equi_sol$w_avg )
Deltaw_fun( equi_sol_no_govt$kappaH, equi_sol_no_govt$kappaL, equi_sol$w_avg )
equi_sol = solve_other_moments(equi_sol_baseline)
equi_sol$private_financing / equi_sol$avg_productivity


# -------------- Table 1 and 2: Calibrated Model Parameters -------------------
desired_order <- c("theta", "lambda_F", "avg_zeta", "lambda",  "beta", "delta", "r", "AH", "AL",
                   "eta_L",  "d_bar", "g_bar")
round_digits <- c(
  theta = 1,
  lambda_F = 1,
  avg_zeta = 2,
  lambda = 2,
  beta = 2,
  delta = 1,
  r = 2,
  AH = 2,
  AL = 2,
  eta_L = 3,
  d_bar = 2,
  g_bar = 2
)
values <- sapply(desired_order, function(name) round(par[[name]], round_digits[name]))
df_parameters <- data.frame(
  Parameter = desired_order,
  Value = values
)
df_moments <- data.frame(
  `Moment Description` = c(
    "Average investment-to-capital ratio",
    "Average GDP drop in crises",
    "Average impact of credit intervention on firm survival likelihood",
    "Average creditor recovery rate",
    "Average output-to-capital ratio",
    "TFP ratio between 90% and 10% percentiles",
    "Average private-sector debt/GDP ratio"
  ),
  Model = c(unlist(calibration_results$key_moments[c(1,2,5,3,4)]), par$AH_AL_ratio, calibration_results$key_moments[[6]]),
  Data = c( moment_targets[c(1,2,5,3,4)], par$AH_AL_ratio, moment_targets[6])
)
write_xlsx(
  list(
    Table1 = df_parameters,
    Table2 = df_moments
  ),
  path = "../Tables/main_tables.xlsx"
)



#-------- Figure 1: Plot the basic solutions --------
q_vec = c(qL, qH)
for(iter in c(1,2)){
  q = q_vec[iter]
  if(q==qL){
    zeta_vec = seq(0,0.3,length.out=800)
  }else{
    zeta_vec = seq(0,0.55,length.out=800)
  }
  
  init_vec = rep(0, length(zeta_vec))
  result = solve_individual_optimal_faster(q, d_bar=d_bar, g_bar)
  result_no_moral_hazard = solve_individual_optimal_faster(q, d_bar=d_bar, g_bar, Moral_Hazard = FALSE)
  result_no_endogenous_limit = solve_individual_optimal_faster(q, d_bar=d_bar, g_bar, endogenous_debt_limit = FALSE)
  result_no_govt = solve_individual_optimal_faster(q, d_bar=d_bar, g_bar=0)
  result_first_best = solve_individual_optimal_faster(q, d_bar=2, g_bar=2, Moral_Hazard=FALSE)
  TB = data.table(zeta=zeta_vec, l=result$l_fun(zeta_vec), l_no_endogenous_limit=result_no_endogenous_limit$l_fun(zeta_vec), 
                  l_first_best=result_first_best$l_fun(zeta_vec),
                  l_no_govt=result_no_govt$l_fun(zeta_vec), x_star=result$endogenous_private_debt_limit_fun(zeta_vec), 
                  continuation=result$C_indicator(zeta_vec), 
                  continuation_no_govt = result_no_govt$C_indicator(zeta_vec),
                  capacity_no_govt = result_no_govt$endogenous_private_debt_limit_fun(zeta_vec),
                  obj_full_payment=result$obj1_fun(zeta_vec), obj_renegotiation = result$obj2_fun(zeta_vec))
  pdf(paste0("../Figures/Figure1_", iter, ".pdf"), width=plot_width, height=plot_height, family="serif")
  par(mar = margins, mfrow=c(1,1))
  MyPlot(TB$zeta, TB[,list(l,l_no_govt, l_first_best)], xlab=TeX("$\\zeta$"), ylab=TeX("financing per unit of capital"), # main=paste0("g_bar=",g_bar,",d_bar=",d_bar),
         LegendNames = c("optimal financing (with intervention)", "optimal financing (laissez-faire)", "first-best financing"),
         LegendPosition = "topright", ylim=c(0,0.5))
  index = first(which(TB$continuation)) # Cutoff for continuation
  index_no_govt =  first(which(TB$continuation_no_govt))
  # abline( v = zeta_vec[index] )
  abline( h=c(d_bar, d_bar+g_bar), lty=3, lwd=1 )
  abline( v= c(zeta_vec[index]) , lty=3, lwd=1) 
  dev.off()
}


# -------------- Figure 2: Government Funding and Normal/Crisis Time Economy ------------
g_bar_vec = seq(0.01, AH , length.out=50)
init_vec = rep(0, length(g_bar_vec))
TB = data.table( g_bar=g_bar_vec, qH=init_vec, qL=init_vec, iH=init_vec, iL=init_vec, creditor_profit_H_frac=init_vec, creditor_profit_L_frac=init_vec, 
                 H_zeta_ = init_vec, L_zeta_ = init_vec, dK=init_vec, dw=init_vec,
                 muw = init_vec, muw_middle=init_vec, muK=init_vec, muK_middle=init_vec, w_avg=init_vec,
                 avg_creditor_recovery = init_vec, avg_investment_to_GDP=init_vec, avg_output_drop=init_vec, takeup_ratio=init_vec,
                 strategic_bankruptcy=init_vec, liquidity_bankruptcy=init_vec, total_bankruptcy=init_vec, 
                 H_zeta_bar = init_vec, L_zeta_bar = init_vec,
                 r1 = init_vec,  r2=init_vec, d_hat1 = init_vec, d_hat2 = init_vec, 
                 xL1 = init_vec, xL2 = init_vec, xH1=init_vec, xH2=init_vec, output_change=init_vec
                 )
zeta_cases = c(0.1, 0.2)   # These two cases map to the interest rate r1 and r2. 
w0=calibration_results$equi_sol$w_steady   # Evaluate at the average w0
# w0 = 0.5  # Evaluate the change of quality and quantity at w=0.5
for(iter in 1:length(g_bar_vec)){
  print( paste0("Progress:" , round(iter/length(g_bar_vec)*100) , "%") )
  equi_sol = list()
  equi_sol = solve_equilibrium(d_bar, g_bar_vec[iter])
  equi_sol = solve_other_moments(equi_sol)
  w0 = equi_sol$w_avg
  TB$qH[iter] = equi_sol$qH
  TB$qL[iter] = equi_sol$qL
  TB$iH[iter] = equi_sol$iH
  TB$iL[iter] = equi_sol$iL
  TB$creditor_profit_H_frac[iter] = equi_sol$creditor_profit_frac_H
  TB$creditor_profit_L_frac[iter] = equi_sol$creditor_profit_frac_L
  TB$H_zeta_[iter] = equi_sol$zeta__H
  TB$L_zeta_[iter] = equi_sol$zeta__L
  TB$dK[iter] = DeltaK_fun( equi_sol$kappaH, equi_sol$kappaL, w0 )
  TB$dw[iter] = Deltaw_fun( equi_sol$kappaH, equi_sol$kappaL, w0 )
  TB$w_avg[iter] = equi_sol$w_avg
  TB$muw[iter] = muw_fun( equi_sol$iH, equi_sol$iL, w0 )
  TB$muw_middle[iter] = muw_fun( equi_sol$iH, equi_sol$iL, 0.5 )
  TB$muK[iter] = muK_fun( equi_sol$iH, equi_sol$iL, w0 )
  TB$muK_middle[iter] = muK_fun( equi_sol$iH, equi_sol$iL, 0.5 )
  TB$avg_creditor_recovery[iter] = equi_sol$avg_creditor_recovery
  TB$avg_investment_to_GDP[iter] = equi_sol$avg_investment_to_GDP
  TB$avg_output_drop[iter] = equi_sol$avg_output_drop
  TB$avg_productivity[iter] = equi_sol$avg_productivity
  TB$takeup_ratio[iter] = equi_sol$takeup_ratio
  TB$strategic_bankruptcy[iter] = equi_sol$strategic_bankruptcy_total
  TB$liquidity_bankruptcy[iter] = equi_sol$liquidity_bankruptcy_total
  TB$total_bankruptcy[iter]     = equi_sol$bankruptcy_total
  TB$H_zeta_bar[iter] = equi_sol$zeta_bar_H
  TB$L_zeta_bar[iter] = equi_sol$zeta_bar_L
  TB$r1[iter] = equi_sol$result_qH$interest_rate_fun(zeta_cases[1])
  TB$r2[iter] = equi_sol$result_qH$interest_rate_fun(zeta_cases[2])
  TB$d_hat1[iter] = equi_sol$result_qL$endogenous_private_debt_limit_fun(zeta_cases[1])
  TB$d_hat2[iter] = equi_sol$result_qL$endogenous_private_debt_limit_fun(zeta_cases[2])
  TB$xL1[iter] = equi_sol$result_qL$l_fun(zeta_cases[1])
  TB$xL2[iter] = equi_sol$result_qL$l_fun(zeta_cases[2])
  TB$xH1[iter] = equi_sol$result_qH$l_fun(zeta_cases[1])
  TB$xH2[iter] = equi_sol$result_qH$l_fun(zeta_cases[2])
  TB$output_change[iter] = Delta_output_fun(equi_sol$kappaH, equi_sol$kappaL, w0)
}
TB_govt_financing = TB

pdf(paste0("../Figures/Figure2a_govt_financing_and_DeltaK.pdf"), width=plot_width, height=plot_height, family="serif")
par(mar = margins, mfrow=c(1,1))
MyPlot(TB$g_bar, TB[,list(dK*100)], xlab=TeX("govt program size $\\bar{g}$"), ylab=TeX("capital quantity change $\\Delta K$  (%)")); abline(h=0)
dev.off()

pdf(paste0("../Figures/Figure2b_govt_financing_and_Deltaw.pdf"), width=plot_width, height=plot_height, family="serif")
par(mar = margins, mfrow=c(1,1))
MyPlot(TB$g_bar, TB$dw*100, xlab=TeX("govt program size $\\bar{g}$"), ylab=TeX("capital quality change  $\\Delta \\omega$ (%)") ); abline(h=0)
dev.off()

TB[,i_diff:=iH-iL]
(max(TB$i_diff) - min(TB$i_diff))/lambda
MyPlot(TB$g_bar, TB[,list(dw)], xlab=TeX("govt program size $\\bar{g}$"), ylab=TeX("$dw$")); abline(h=0)

#---------------- Figure 3: Government Lending and Private Lending -------------
pdf(paste0("../Figures/Figure3a_govt_financing_and_xH_and_xL.pdf"), width=plot_width, height=plot_height, family="serif")
par(mar = margins, mfrow=c(1,1))
MyPlot(TB$g_bar, TB[,list(xL1,xH1)],  xlab=TeX("govt program size $\\bar{g}$"), ylab=TeX("optimal financing"), LegendNames = c("L type (strategic default)", "H type") , abline_h=0 , ylim=c(0,0.7))
dev.off()

pdf(paste0("../Figures/Figure3b_govt_financing_and_decomposing_xL.pdf"), width=plot_width, height=plot_height, family="serif")
par(mar = margins, mfrow=c(1,1))
MyPlot(TB$g_bar, TB[,list(xL1, d_hat1, g_bar)],  xlab=TeX("govt program size $\\bar{g}$"), ylab=TeX("optimal financing of L type"), 
       LegendNames = c("total borrowing", TeX("private-sector financing"), TeX("government financing") ) ,  abline_h=0 )
dev.off()


# -------------------- Figure 4: Welfare -------------------------
# Solve for the Welfare function
g_bar_vec = seq(0,AH,length.out=30)
w_grid = seq(0.0001,0.9999,length.out = 500)
init_vec = rep(0, length(g_bar_vec))
index = length(w_grid)/2+1
# w0 = w_grid[index]
TB = data.table( g_bar=g_bar_vec, g_takeup_avg=init_vec, g_takeup=init_vec, g_takeup2=init_vec, g_takeup3=init_vec, productivity_avg=init_vec, welfare_avg=init_vec, welfare = init_vec, welfare2 = init_vec, welfare3=init_vec, l_dN=init_vec, w0=init_vec, muw=init_vec, muK=init_vec, Delta_w=init_vec, Delta_K=init_vec, qH=init_vec, qL=init_vec )

for(iter in 1:length(g_bar_vec)){
  print( paste0("Progress: ", round(100*iter/length(g_bar_vec),1), "%") )
  result_welfare = solve_welfare(w_grid, d_bar, g_bar=g_bar_vec[iter])
  # Another idea: compare the welfare with and without the slippery slope of credit intervention. 
  equi_sol = solve_other_moments(result_welfare$equi_sol)
  w0 = equi_sol$w_avg
  TB$welfare[iter] = result_welfare$W_fun(0.01)
  TB$welfare2[iter] = result_welfare$W_fun(0.5)
  TB$welfare3[iter] = result_welfare$W_fun(1)
  TB$welfare_avg[iter] = result_welfare$W_fun(w0)
  TB$g_takeup[iter] = equi_sol$total_take_up_fun(0.01)
  TB$g_takeup2[iter] = equi_sol$total_take_up_fun(0.5)
  TB$g_takeup3[iter] = equi_sol$total_take_up_fun(0.99)
  TB$g_takeup_avg[iter] = equi_sol$total_takeup
  TB$productivity_avg[iter] = equi_sol$avg_productivity
  TB$l_dN[iter] = interp1( result_welfare$w_grid, result_welfare$l_dN, w0 )
  TB$w0[iter] = w0
  TB$muw[iter]  = interp1( result_welfare$w_grid, result_welfare$muw, w0 )
  TB$muK[iter] = interp1( result_welfare$w_grid, result_welfare$muK, w0 ) 
  TB$Delta_w[iter] = interp1( result_welfare$w_grid, result_welfare$Delta_w, w0 )
  TB$Delta_K[iter] = interp1( result_welfare$w_grid, result_welfare$Delta_K, w0 ) 
  TB$qH[iter] = result_welfare$qH
  TB$qL[iter] = result_welfare$qL
}
pdf(paste0("../Figures/Figure4a_w_avg_and_g_bar.pdf"), width=plot_width, height=plot_height, family="serif")
par(mfrow=c(1,1), mar=margins)
MyPlot(TB$g_bar, TB[,list(w0)], xlab=TeX("govt program size $\\bar{g}$"), ylab=TeX("average capital quality") )
dev.off()

index = which(TB$welfare_avg==max(TB$welfare_avg))
pdf(paste0("../Figures/Figure4b_welfare_and_g_bar.pdf"), width=plot_width, height=plot_height, family="serif")
par(mfrow=c(1,1), mar=margins)
MyPlot(TB$g_bar, TB[,list(100*(welfare_avg/welfare_avg[1]-1))], xlab=TeX("govt program size $\\bar{g}$"), ylab=TeX("welfare (% diff from laissez-faire)"), abline_v=TB$g_bar[index], abline_h=0  )
dev.off()

# ----- Figure 5: Optimal Static Welfare ----------
g_vec = seq(0.15,0.35,length.out=150)
w_grid = w_grid=seq(0,1,length.out=50)
plot_process = TRUE
result_optimal_W = Optimal_static_welfare(calibration_results, w_grid=w_grid, g_vec = g_vec , plot_process = TRUE)

TB = data.table(w=w_grid,g_bar_optimal=result_optimal_W$g_bar_optimal, govt_lending_takeup=rep(0,length(w_grid)),  W_optimal=result_optimal_W$W_optimal, W_optimal_shutoff_slippery_slope=result_optimal_W$W_optimal_shutoff_slippery_slope,
                g_optimal_shutoff_slippery_slope=result_optimal_W$g_optimal_shutoff_slippery_slope_vec, govt_lending_takeup_shutoff_slippery_slope=rep(0,length(w_grid))  )
for(iter in 1:nrow(TB)){
  g_bar_optimal = TB$g_bar_optimal[iter]
  g_bar_optimal_shutoff_slippery_slope = TB$g_optimal_shutoff_slippery_slope[iter]
  equi_sol = calibration_results$f_list$solve_equilibrium(d_bar=calibration_results$equi_sol$d_bar,g_bar=g_bar_optimal)
  equi_sol = calibration_results$f_list$solve_other_moments(equi_sol)
  equi_sol_shutoff_slippery_slope = calibration_results$f_list$solve_equilibrium(d_bar=calibration_results$equi_sol$d_bar,g_bar=g_bar_optimal_shutoff_slippery_slope)
  equi_sol_shutoff_slippery_slope = calibration_results$f_list$solve_other_moments(equi_sol_shutoff_slippery_slope)
  TB$govt_lending_takeup[iter] = equi_sol$total_takeup
  TB$govt_lending_takeup_shutoff_slippery_slope[iter] = equi_sol_shutoff_slippery_slope$total_takeup
}
# Next, plot the optimal solution without
pdf(paste0("../Figures/Figure5a_optimal_g_bar.pdf"), width=plot_width, height=plot_height, family="serif")
par(mfrow=c(1,1), mar=margins)
MyPlot(TB$w, TB[,list(g_bar_optimal, g_optimal_shutoff_slippery_slope)], 
       LegendNames = c(TeX("Optimal $\\bar{g}^*$"), TeX("Pseudo-optimal $\\tilde{g}$")), 
       xlab=TeX("$\\omega_0$"), ylab=TeX("govt program size $\\bar{g}$") )
dev.off()
# Welfare gap between the economy realizing the slippery slope and the economy that does not. 
pdf(paste0("../Figures/Figure5b_optimal_g_bar_welfare_gap.pdf"), width=plot_width, height=plot_height, family="serif")
par(mfrow=c(1,1), mar=margins)
MyPlot(TB$w, TB[,list(100*(W_optimal_shutoff_slippery_slope-W_optimal)/W_optimal)],  
       xlab=TeX("$\\omega_0$"), ylab="welfare gap (%)", abline_h = 0)
dev.off()

# -------------------- Figure 6: Slippery Slope in the Same Economy --------------------
# Government intervention and the avg w
g_bar_vec = c( seq(0, g_bar ,length.out=20), seq(g_bar+0.005, 0.25, length.out=20) )
index_baseline = which.min( abs(g_bar_vec-0.14)  )
GDP_drop_target = -0.1
# intervention and next intervention scale. 
govt_funding_next_crisis_vec = rep(0, length(g_bar_vec))
equi_sol_no_govt = solve_equilibrium(d_bar, g_bar=0)
equi_sol_no_govt = solve_other_moments(equi_sol_no_govt)
w0 = equi_sol_no_govt$w_avg
K0 = 1
# Think about the actual amount of intervention. 
govt_credit_provision_vec = g_bar_vec
govt_credit_provision_next_crisis = g_bar_vec
w_next_crisis_vec = g_bar_vec
K_next_crisis_vec = g_bar_vec
productivity_vec = rep(0,length(g_bar_vec))
productivity_next_crisis_vec = rep(0,length(g_bar_vec))
for(iter_g_bar in 1:length(g_bar_vec)){
  print( paste0("Progress:", iter_g_bar/length(g_bar_vec)*100, "%") )
  resultH = post_crisis_fraction_of_capital(equi_sol_no_govt$qH, d_bar, g_bar_vec[iter_g_bar], type="H")
  resultL = post_crisis_fraction_of_capital(equi_sol_no_govt$qL, d_bar, g_bar_vec[iter_g_bar], type="L")
  kappaH = resultH$kappa
  kappaL = resultL$kappa
  govt_credit_provision_vec[iter_g_bar] = w0*resultH$govt_sector_financing + (1-w0)*resultL$govt_sector_financing
  productivity_vec[iter_g_bar] = w0*AH + (1-w0)*AL
  w = w0  + Deltaw_fun(kappaH, kappaL, w0) # immediately after the crisis
  K = K0 * (1 + DeltaK_fun(kappaH, kappaL, w0) )
  for(iter in 1:10/dt){
    # 10 years
    w = w + muw_fun(equi_sol_no_govt$iH, equi_sol_no_govt$iL, w)*dt
    K = K * (1 + muK_fun(equi_sol_no_govt$iH, equi_sol_no_govt$iL, w)*dt )
  }  
  w_next_crisis_vec[iter_g_bar] = w
  K_next_crisis_vec[iter_g_bar] = K
  result = uniroot( function(g_bar) output_drop_and_gbar(g_bar, w, equi_sol_no_govt$qH, equi_sol_no_govt$qL) - GDP_drop_target, interval = c(0, 2)  )
  govt_funding_next_crisis_vec[iter_g_bar] = result$root
  
  resultH= post_crisis_fraction_of_capital(equi_sol_no_govt$qH, d_bar, govt_funding_next_crisis_vec[iter_g_bar], type="H")
  resultL= post_crisis_fraction_of_capital(equi_sol_no_govt$qL, d_bar, govt_funding_next_crisis_vec[iter_g_bar], type="L")
  # Total credit provision
  govt_credit_provision_next_crisis[iter_g_bar] = w*resultH$govt_sector_financing + (1-w)*resultL$govt_sector_financing
  productivity_next_crisis_vec[iter_g_bar] = w*AH + (1-w)*AL
}
sensitivity = ( govt_funding_next_crisis_vec[index_baseline] - govt_funding_next_crisis_vec[1] ) / ( g_bar_vec[index_baseline] - g_bar_vec[1] )
pdf(paste0("../Figures/Figure6a_same_economy_slippery_slope_w.pdf"), width=plot_width,height=plot_height, family="serif")
par(mfrow=c(1,1), mar=margins)
MyPlot( g_bar_vec, w_next_crisis_vec, xlab=TeX("current govt program size $\\bar{g}$"), ylab=TeX("next crisis (in 10 years) $\\omega$'") )
dev.off()
pdf(paste0("../Figures/Figure6b_same_economy_slippery_slope_credit.pdf"), width=plot_width,height=plot_height, family = "serif")
par(mfrow=c(1,1), mar=margins)
MyPlot( g_bar_vec, govt_funding_next_crisis_vec, xlab=TeX("current govt program size $\\bar{g}$"), ylab=TeX("next crisis (in 10 years) govt intervention $\\bar{g}'$") )
dev.off()

# ----- Figure 7: Slippery Slope Across Two Economies ----------
w_starting_vec = equi_sol_no_govt$w_steady 
years_sim = 20
GDP_drop_target = -0.1
# Set expectation as non-government intervention!!! 
for(w_starting in w_starting_vec){
  dt = 1/4
  # dN_vec = c( rep(0,10/dt), 1, rep(0,(years_sim-10)/dt-1) )
  dN_vec = c( 1, rep(0,10/dt), 1, rep(0,(years_sim-10)/dt-2) )
  # dN_vec = rep( 0, years_sim/dt )
  N = first(which(dN_vec==1))
  equi_sol_no_govt = solve_equilibrium( d_bar, g_bar=0 )
  # par1: the parameters for case 1, i.e., expectation is with intervention and there is actually government intervention. The intervention gamma is given by low gamma (lenient intervention)
  equi_sol_with_govt = solve_equilibrium( d_bar, g_bar )
  equi_sol_with_govt$qH = equi_sol_no_govt$qH
  equi_sol_with_govt$qL = equi_sol_no_govt$qL
  equi_sol_with_govt$iH = equi_sol_no_govt$iH
  equi_sol_with_govt$iL = equi_sol_no_govt$iL
  
  simulation_case1 = simulation_fun( equi_sol_with_govt, dN_vec = dN_vec  , dt=dt , years_sim = years_sim, w_init=w_starting)
  simulation_case1[  , actual_govt_intervention := 0 ]  # scale of govt intervention
  simulation_case1[ which(dN_vec==1) , actual_govt_intervention := equi_sol_with_govt$total_borrowing_H*w + equi_sol_with_govt$total_borrowing_L*(1-w)  ]  # scale of govt intervention
  simulation_case1[  , gbar_needed := 0  ]
  simulation_case1[  , total_govt_financing := 0  ]
  simulation_case1[  , total_private_financing := 0  ]
  # par2: the parameters for case 2, i.e., expectation is laissez-faire, and there is laissez-faire at first crisis, but with intervention at the second crisis to achieve the same GDP decline
  simulation_case2 = simulation_fun( equi_sol_no_govt, dN_vec = dN_vec, dt=dt , years_sim = years_sim, w_init=w_starting)
  # Then solve for the gamma that achieves the same GDP decline.
  simulation_case2[  , gbar_needed := 0  ]
  simulation_case2[  , total_govt_financing := 0  ]
  simulation_case2[  , total_private_financing := 0  ]
  
  MyPlot( zeta_vec, equi_sol_with_govt$result_qL$l_fun(zeta_vec) )
  
  for(iter in 1:nrow(simulation_case2)){
    print( paste0( "Progress:", round( iter/nrow(simulation_case2)*100, 0 )  , "%" ) )
    w = simulation_case1[iter]$w
    result = uniroot( function(g_bar) output_drop_and_gbar(g_bar, w, equi_sol_no_govt$qH, equi_sol_no_govt$qL) - GDP_drop_target, interval = c(0, 2)  )
    g_bar_needed = result$root
    simulation_case1[iter]$gbar_needed = g_bar_needed
    resultH = post_crisis_fraction_of_capital( equi_sol_no_govt$qH, equi_sol_no_govt$d_bar, g_bar_needed, type="H"  )  # Then find the total governmetn credit provision
    resultL = post_crisis_fraction_of_capital( equi_sol_no_govt$qL, equi_sol_no_govt$d_bar, g_bar_needed, type="L"  )
    simulation_case1[iter]$total_govt_financing = w*resultH$govt_sector_financing + (1-w)*resultL$govt_sector_financing
    simulation_case1[iter]$total_private_financing = w*resultH$private_sector_financing + (1-w)*resultL$private_sector_financing
    
    w = simulation_case2[iter]$w
    result = uniroot( function(g_bar) output_drop_and_gbar(g_bar, w, equi_sol_no_govt$qH, equi_sol_no_govt$qL) - GDP_drop_target, interval = c(0, 2), tol=10^(-12)  )
    g_bar_needed = result$root
    simulation_case2[iter]$gbar_needed = g_bar_needed
    resultH = post_crisis_fraction_of_capital( equi_sol_no_govt$qH, equi_sol_no_govt$d_bar, g_bar_needed, type="H"  )  # Then find the total governmetn credit provision
    resultL = post_crisis_fraction_of_capital( equi_sol_no_govt$qL, equi_sol_no_govt$d_bar, g_bar_needed, type="L"  )
    simulation_case2[iter]$total_govt_financing = w*resultH$govt_sector_financing + (1-w)*resultL$govt_sector_financing
    simulation_case2[iter]$total_private_financing = w*resultH$private_sector_financing + (1-w)*resultL$private_sector_financing
  }
  pdf(paste0("../Figures/Figure7a_two_economy_slippery_slope_path_of_w.pdf"), width=plot_width, height=plot_height, family="serif")
  par(mfrow=c(1,1), mar=margins)
  MyPlot( simulation_case1$T, data.frame( simulation_case1$w, simulation_case2$w ), xlab="simulation years", ylab=TeX("firm quality $\\omega_t$"), LegendNames = c("with intervention", "laissez-faire") )
  dev.off()
  
  pdf(paste0("../Figures/Figure7b_two_economy_slippery_slope.pdf"), width=plot_width, height=plot_height, family="serif")
  par(mfrow=c(1,1), mar=margins)
  MyPlot( simulation_case1$T, 100*(simulation_case1$gbar_needed-simulation_case2$gbar_needed )/simulation_case2$gbar_needed, xlab="simulation years", ylab="extra govt intervention needed (%)" )
  dev.off()
}




# ----- Figure A1: Optimal Dynamic Policy ----------
# Note: result_optimal_W is the optimal welfare in the static govt intervention case. It is used below to initialize the algorithm for optimal dynamic welfare case. 
result_dynamic_optimal_welfare = solve_optimal_welfare_dynamic_gbar(calibration_results, result_optimal_W)
# Q: is the dynamic welfare algorithm converging well?
MyPlot(1:length(result_dynamic_optimal_welfare$diff_vec), log(result_dynamic_optimal_welfare$diff_vec))
# This contains the equilibrium solution
result_dynamic_gbar_equi = result_dynamic_optimal_welfare$result_dynamic_equi   
# Store the grid and optimal welfare. 
w_grid = result_dynamic_optimal_welfare$w_grid
W_dynamic_optimal = result_dynamic_optimal_welfare$result_welfare$W_fun(w_grid)
W_static_optimal = approx(result_optimal_W$w_grid,result_optimal_W$W_optimal,w_grid)$y

# Then plot the welfare gap between the dynamic optimizing case and the static optimizing case.
pdf(paste0("../Figures/FigureA1a_optimal_dynamic_g_bar.pdf"), width=plot_width, height=plot_height, family="serif")
par(mfrow=c(1,1), mar=margins)
MyPlot( w_grid,  data.frame(approx(result_optimal_W$w_grid,result_optimal_W$g_bar_optimal,w_grid)$y, result_dynamic_optimal_welfare$g_vec_optimal ), 
        LegendNames = c(TeX("static optimal"),TeX("dynamic optimal")), xlab=TeX("$\\omega$"), ylab=TeX("$\\bar{g}$") )
dev.off()
pdf(paste0("../Figures/FigureA1b_optimal_dynamic_W.pdf"), width=plot_width, height=plot_height, family="serif")
par(mfrow=c(1,1), mar=margins)
MyPlot( w_grid,  100*data.frame(W_dynamic_optimal/W_static_optimal-1 ), ylab = "welfare difference (%)",  xlab=TeX("$\\omega$") )
dev.off()
# Suppose that we use the dynamic g_bar policy that maximizes a static case. Will welfare be the same at w=0?
result_dynamic_comparison = solve_for_equilibrium_with_dynamic_g(calibration_results, g_fun = approxfun(result_optimal_W$w_grid, result_optimal_W$g_bar_optimal) )
W_dynamic_comparison = solve_welfare_dynamic_gbar(result_dynamic_comparison, calibration_results)
par(mfrow=c(1,1), mar=margins)
MyPlot( w_grid,  100*data.frame(W_dynamic_comparison$W_fun(result_optimal_W$w_grid)/W_static_optimal-1 ), ylab = "Welfare difference (%)",  xlab=TeX("$\\omega$") )
MyPlot( w_grid,  100*data.frame(result_optimal_W$g_bar_optimal/result_optimal_W$g_bar_optimal-1 ), ylab = "Welfare difference (%)",  xlab=TeX("$\\omega$") )
TB = as.data.table(result_dynamic_gbar_equi$TB)

# ----- Figure A2: Slippery Slope Under Optimal Dynamic Policy ----------
dt = 1/4
years_sim = 20
equi_sol_no_govt = solve_other_moments(equi_sol_no_govt)
w_starting = equi_sol_no_govt$w_steady
dN_vec = c( 1, rep(0,10/dt), 1, rep(0,(years_sim-10)/dt-2) )
# Dynamic optimal intervention. 
simulation_case1 = simulation_dynamic_gbar_fun( result_dynamic_gbar_equi, dN_vec = dN_vec  , dt=dt , years_sim = years_sim, w_init=w_starting)
# Note: simulation results include the GDP drop in the next immediate crisis, GDP_drop_in_crisis
# Static optimal intervention: laissez-faire and no such expectation.
equi_sol_no_govt = solve_equilibrium( d_bar, g_bar=0 )
simulation_case2 = simulation_fun( equi_sol_no_govt, dN_vec = dN_vec, dt=dt , years_sim = years_sim, w_init=w_starting)
# Then solve for the gamma that achieves the same GDP decline.
simulation_case2[  , g_bar := 0  ]
simulation_case2[  , govt_credit := 0  ]
# Then solve for the required intervention in the laissez-faire case. 
for(iter in 1:nrow(simulation_case2)){
  w = simulation_case2$w[iter]
  print( paste0( "Progress:", round( iter/nrow(simulation_case2)*100, 0 )  , "%" ) )
  result = uniroot( function(g_bar) output_drop_and_gbar(g_bar, w, equi_sol_no_govt$qH, equi_sol_no_govt$qL) - simulation_case1$GDP_drop_in_crisis[iter], interval = c(0, 1)  )
  gbar_needed = result$root
  simulation_case2[iter]$g_bar = gbar_needed
  result_H = post_crisis_fraction_of_capital(equi_sol_no_govt$qH, par$d_bar,  g_bar = gbar_needed, type="H")
  result_L = post_crisis_fraction_of_capital(equi_sol_no_govt$qL, par$d_bar, g_bar = gbar_needed, type="L")
  simulation_case2$govt_credit[iter] = w*result_H$govt_sector_financing + (1-w)*result_L$govt_sector_financing 
}
pdf(paste0("../Figures/FigureA2b_two_economy_dynamic_optimal_slippery_slope.pdf"), width=plot_width, height=plot_height, family="serif")
par(mfrow=c(1,1), mar=margins)
MyPlot( simulation_case1$T, 100*(simulation_case1$g_bar-simulation_case2$g_bar )/simulation_case2$g_bar , xlab="simulation years", ylab="extra govt financing needed (%)")
dev.off()
pdf(paste0("../Figures/FigureA2a_two_economy_dynamic_optimal_slippery_slope_path_of_w.pdf"), width=plot_width, height=plot_height, family="serif")
par(mfrow=c(1,1), mar=margins)
MyPlot( simulation_case1$T, data.frame( simulation_case1$w, simulation_case2$w ), xlab="simulation years", ylab=TeX("firm quality $\\omega_t$"), LegendNames = c("optimal dynamic govt intervention", "laissez-faire") )
dev.off()

# ------- Figure A3: Govt financing and capital values ---------
pdf(paste0("../Figures/FigureA3a_govt_financing_and_qL.pdf"), width=plot_width, height=plot_height, family="serif")
par(mar = margins, mfrow=c(1,1))
MyPlot(TB_govt_financing$g_bar, TB_govt_financing[,list(qL)], xlab=TeX("govt program size $\\bar{g}$"), ylab=TeX("$q^L$")); abline(h=0)
dev.off()

pdf(paste0("../Figures/FigureA3b_govt_financing_and_qH.pdf"), width=plot_width, height=plot_height, family="serif")
par(mar = margins, mfrow=c(1,1))
MyPlot(TB_govt_financing$g_bar, TB_govt_financing[,list(qH)], xlab=TeX("govt program size $\\bar{g}$"), ylab=TeX("$q^H$")); abline(h=0)
dev.off()

# ----- Figure A4: Stationary Distribution ----------
sim_N = 10000/dt # 1000 years of simulation
g_bar_vec = c(0,g_bar)
equi_sol_list = list()
TBs = list()
for(iter_gbar in 1:length(g_bar_vec)){
  equi_sol = solve_equilibrium(d_bar, g_bar_vec[iter_gbar])
  equi_sol_list[[iter_gbar]] = equi_sol
  TB = Simulate_model(sim_N, seed=1, d_bar, g_bar=g_bar_vec[iter_gbar])
  par(mar = margins, mfrow=c(1,2))
  index = 1:(200/dt)
  MyPlot(TB[index]$t, TB[index]$w, xlab="years", ylab=TeX("\\omega"), main=TeX( paste0("govt program size $\\bar{g}=$", g_bar_vec[iter_gbar]))  )
  MyPlot(TB[index]$t, log(TB[index]$K), xlab="years", ylab="log(K)", main=paste0("g_bar=", g_bar_vec[iter_gbar]))
  TBs[[iter_gbar]] = TB
}
TBs[[1]][,category:="laissez-faire"]
TBs[[2]][,category:="baseline intervention"]
TB = rbind( TBs[[1]], TBs[[2]] )
TB[,category:=as.factor(category)]
pdf(paste0("../Figures/FigureA4_stationary_distribution.pdf"), width=8, height=4, family = "serif")
ggplot(TB, aes(x=w, fill=category, color=category)) +  geom_density(alpha=0.3, bw=0.008) + xlab(TeX("capital quality $\\omega$")) + theme(plot.title = element_text(hjust = 0.5))
dev.off()




